% This file calculates the population shocks between 1950 and 2015
%

clear

age_profile_00_to_49 = csvread('../data/demographic/age_profile_00_to_49.csv',1,1) ;
age_profile_all      = csvread('../data/demographic/age_profile.csv',1,1) ;
age_profile_all      = [age_profile_00_to_49 age_profile_all] ;

age = (0:85)' ;
idxa = 1:length(age) ;

timevec = 1900:2015 ;
idxt = 1:length(timevec) ;

age_profile = age_profile_all( ...
    idxa(age==16):idxa(age==85), ...
    idxt(timevec==1950):idxt(timevec==2015)) ;

clear gam_mat

gam_probs_all = csvread('../data/demographic/actuarial_table_over_time.csv',1,1) ;

age = (0:119)' ;
idxa = 1:length(age) ;

timevec = 1900:10:2070 ;
idxt = 1:length(timevec) ;

timevec_new = 1860:2070 ;
idxt_new = 1:length(timevec_new) ;

for jj=1:length(age)
    gam_mat(jj,:) = ...
        interp1(timevec', gam_probs_all(jj,:)', timevec_new','linear','extrap') ;
end

gam_mat(gam_mat>1) = 1 ;

gam_sh_mat = gam_mat(idxa(age==16):idxa(age==95),find(timevec_new==1940):find(timevec_new==2070)) - ...
    repmat(gam_mat(idxa(age==16):idxa(age==95),find(timevec_new==2070)),[1 (2070-1940+1)]) ;

gam_mat_ = zeros(size(gam_mat));
for i=1:length(gam_mat(:,1))
    gam_mat_(i,:) = prod(1-gam_mat(1:i,:));
end
life_exp = sum(gam_mat_);

gam_mat_b = zeros(size(gam_mat));
for i=65:length(gam_mat(:,1))
    gam_mat_b(i,:) = prod(1-gam_mat(65:i,:));
end

life_exp_at_65 = sum(gam_mat_b);

gam_mat_orig   = gam_mat ;

gam_mat = gam_mat(idxa(age==16):idxa(age==95),find(timevec_new==1940):find(timevec_new==2070)) ;

clear gam_mat_new

for ii=1:length(gam_mat(:,1))
	gam_mat_new(ii,:) = ...
        gam_mat_orig(idxa(age==(16+ii-1)),(find(timevec_new==1940)-ii+1):(find(timevec_new==1940)-ii+length(gam_mat(1,:)))) ;
end

gam_mat = gam_mat_new ;

syn_n = (age_profile(1,:)/age_profile(1,1))';

% want syn_n starts in 1949
for ii=2:length(syn_n)
    syn_n(ii) = syn_n(ii-1)*(1-gam_mat(ii-1,10)) ;
end

n0_data_sh = age_profile(1,:)./sum(age_profile) ;

nsh(1) = 1 ;

for tt=2:length(age_profile(1,:))
    for jj=2:length(age_profile(:,1))
        syn_n(jj,tt) = syn_n(jj-1,tt-1)*(1-gam_mat(jj-1,tt-1+10)) ;
    end
    
    nsh(tt) = n0_data_sh(tt)*(sum(syn_n(2:end,tt))) / (1-n0_data_sh(tt)) ;
    
    syn_n(1,tt) = nsh(tt) ;
end

for tt=length(age_profile(1,:))+1:(131-7)
    nsh(tt) = (nsh(tt-1)-1)*0.9 + 1 ;
end

nsh = nsh-1 ;

nshocks = nsh(1:60)' ;

% save shocks here
%save nsh nshocks
